{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survival Analysis\n",
    "\n",
    "This notebook presents code and exercises from Think Bayes, second edition.\n",
    "\n",
    "Copyright 2016 Allen B. Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "% matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "from thinkbayes2 import Pmf, Cdf, Suite, Joint\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Weibull distribution\n",
    "\n",
    "The Weibull distribution is often used in survival analysis because it models the distribution of lifetimes for manufactured products, at least over some parts of the range.\n",
    "\n",
    "The following functions evaluate its PDF and CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def EvalWeibullPdf(x, lam, k):\n",
    "    \"\"\"Computes the Weibull PDF.\n",
    "\n",
    "    x: value\n",
    "    lam: parameter lambda in events per unit time\n",
    "    k: parameter\n",
    "\n",
    "    returns: float probability density\n",
    "    \"\"\"\n",
    "    arg = (x / lam)\n",
    "    return k / lam * arg**(k-1) * np.exp(-arg**k)\n",
    "\n",
    "def EvalWeibullCdf(x, lam, k):\n",
    "    \"\"\"Evaluates CDF of the Weibull distribution.\"\"\"\n",
    "    arg = (x / lam)\n",
    "    return 1 - np.exp(-arg**k)\n",
    "\n",
    "def MakeWeibullPmf(lam, k, high, n=200):\n",
    "    \"\"\"Makes a PMF discrete approx to a Weibull distribution.\n",
    "\n",
    "    lam: parameter lambda in events per unit time\n",
    "    k: parameter\n",
    "    high: upper bound\n",
    "    n: number of values in the Pmf\n",
    "\n",
    "    returns: normalized Pmf\n",
    "    \"\"\"\n",
    "    xs = np.linspace(0, high, n)\n",
    "    ps = EvalWeibullPdf(xs, lam, k)\n",
    "    return Pmf(dict(zip(xs, ps)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy also provides functions to evaluate the Weibull distribution, which I'll use to check my implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33093633846922332"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import weibull_min\n",
    "\n",
    "lam = 2\n",
    "k = 1.5\n",
    "x = 0.5\n",
    "\n",
    "weibull_min.pdf(x, k, scale=lam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33093633846922332"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EvalWeibullPdf(x, lam, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1175030974154046"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weibull_min.cdf(x, k, scale=lam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.11750309741540454"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EvalWeibullCdf(x, lam, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what the PDF looks like, for these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEPCAYAAACHuClZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfWd//HXJ/tCAmHfww6CsqksbkQRBTcc27pW60yn\ntdPaOnYWu8yM2mn70860Hf0502prF/vr1Lq0ShUBFSICgiCy76BhCZtsCSQh2/f3x705uYlZSW7O\nzT3v5+NxH5zvueec+0nE++G7m3MOERGRtkjwOwAREen8lExERKTNlExERKTNlExERKTNlExERKTN\nlExERKTNop5MzGy2mW0zsx1m9lAj1zxpZjvNbJ2ZTQyfG2hmi81ss5ltNLNvRFyfY2aLzGy7mS00\ns67R/jlERKRxUU0mZpYAPAVcC4wD7jCzMfWumQMMd86NBO4Dfh5+qxL4pnNuHDAd+FrEvd8C3nLO\njQYWA9+O5s8hIiJNi3bNZAqw0zlX4JyrAJ4H5ta7Zi7wHIBzbhXQ1cz6OOcOOefWhc+fBrYCAyLu\n+W34+LfAzdH9MUREpCnRTiYDgH0R5f3UJoTGrjlQ/xozGwJMBFaGT/V2zh0GcM4dAnq3W8QiItJq\nMd8Bb2ZdgJeAB5xzZxq5TGvCiIj4KCnKzz8ADI4oDwyfq3/NoIauMbMkQonkd865VyOuORxuCjts\nZn2BIw19uJkpyYiInAPnnLXm+mjXTFYDI8ws18xSgNuBefWumQfcA2Bm04CTNU1YwK+ALc65Jxq4\n597w8ReAV2mEc04v53j44Yd9jyFWXvpd6Heh30XTr3MR1ZqJc67KzO4HFhFKXM8657aa2X2ht90z\nzrn5Znadme0CzhBOEmZ2KXAXsNHMPiTUlPUd59wC4HHgBTP7G6AAuDWaP4eIiDQt2s1chL/8R9c7\n93S98v0N3LccSGzkmceBq9sxTBERaYOY74CX9pGXl+d3CDFDv4ta+l3U0u+ibexc28c6AzNz8fzz\niYhEg5nhYqwDXkREAkDJRERE2kzJRERE2kzJRERE2kzJRERE2kzJRERE2izqkxaD5EzpWd56bxub\ndh5g/6ETVDtHTnYGY4b25ZJJwxmZ2xuzVo22ExHpFDTPpB1UV1fz57fX8fKiDzlbXtHodWOH9+Oe\nudMYmdsn6jGJiJyrc5lnomTSRsVnyvjJb95iw479LYsJ+Mw1k7ltzkUkJKiVUURij5JJPdFOJiWl\n5Tz8339hz76j3rlBfXO47ooLGDOsL8lJiRQeOcmKdXtYumYn1dXV3nUXjBrAQ1+8lvS0lKjFJyJy\nLpRM6olmMqmsrOJ7P3udzbsKQ58F3DIrVONITPx0jePg0VM8/cJSNu6o3c5l+KBe/OvfXU9WZlpU\nYhQRORdKJvVEM5n8bt5KXnl7nVf+ym1XMOuSsU3eU11dzQsLP+DFBR9454YN6sX37r9RNRQRiRla\nm6uDbN5VyKsRieSO66c0m0gAEhISuH3Oxdx36xXU/Ffas+8oP3p2EZWVVVGKVkQk+pRMWqm8opKn\nfr/E23R+/KiBfGbWpFY945pLx3LfbVd45Q079vOrP61oxyhFRDqWkkkrzV+6iSPHiwHITE/l/rvy\nzmnuyKxLxnLbnIu88sLlm1m8clu7xSki0pGUTFqh+EwZLy9a65Vvv+4ienTrcs7P+9y1F3Lp5BFe\n+ekX36Wg8FibYhQR8YOSSSv86c0PKSkrB6Bfr65c04J+kqaYGV+9fQaD+3UHQiPEfvrc25RXVLY5\nVhGRjqRk0kJnSs+ycPkWr3zXDVNJSmpwi/pWSUtN5pv3ziI5/Kx9B4/z+7+83+bnioh0JCWTFnpz\nxVZvqZSBfXKYNmFouz17UN8c7r35Eq/8+jsb2FVwpN2eLyISbUomLVBZWcX8pRu98k1XjW/3BRuv\nvWwsE0YPBMABP/vjUqqqqpu+SUQkRiiZtMCqjR9z7OQZALK7pHP5hSPb/TPMjC997nKvuevjA5/w\nekQCExGJZUomLbBkVe2Q3WsvG0tKcnRW7u/Xqyufm32hV/7D66u9YcgiIrFMyaQZx0+dYd3WfV75\nqqljovp5c6+cwKC+OUBoguQvX1xGPC95IyLxQcmkGUvX7PRmu48b0Z/e3bOi+nlJSYl85bYZXvmD\nLQWs3bI3qp8pItJWSiZNcM6xZNV2r3zllNEd8rljhvVl5rTaGtDv5q1UZ7yIxDQlkybsO3SC/YdP\nAJCSnMT0icM67LPvuH4KqSnJXhxL3t/ezB0iIv5RMmnCqg0feccXjsslLTW5wz47JzuDm2dO8MrP\nz19N2dnGtwQWEfGTkkkT3t/4sXc8bXz7TVJsqZuunEBOdgYAJ4pKeHXx+g6PQUSkJZRMGnH0eLG3\nHW9iYgKTxg7q8BjSUpO54/qLvfKri9dzoqikw+MQEWmOkkkjImslF4wcQGZ6qi9xXDlltLcQ5Nny\nijq7NIqIxAolk0as2VTgHU+5YIhvcSQkJHD3TdO88lsrt3JUExlFJMYomTSgvKKSLXsOeuULx+X6\nGA1MOm8Qo4f2BaCqqpqX31zbzB0iIh1LyaQBW/cc8vZkH9C7Gz1zzn0DrPZgZnV2ZVy8aruWWRGR\nmKJk0oD122qXT5kwZqCPkdQaP2oAY4ZF1E4WqXYiIrFDyaQB67cf8I4vGBUbycTMuH1O7ciuxau2\nc/hYkY8RiYjUUjKp51RxKR8f+ASABDPOH9Hf54hqnT+yP2OH9wOgurqalxaqdiIisUHJpJ5Nuwq9\n45FD+pCRnuJjNHWZGbdfV1s7yV+9QyO7RCQmKJnUs3V37SiuC0YN8DGSho0bUbd2Mm+JZsWLiP+U\nTOrZuueQd3xeuMM71twya7J3/OaKrZwqLvUxGhERJZM6SsvKKQj3lxgwKrePvwE1YuKYgQwZ0BOA\nisoqXn9H2/uKiL+UTCLsKDjibYQ1uH+PmOoviWRm3DJrkld+491NlJSW+xiRiASdkkmEbRFNXGOG\nxmYTV43pE4bSv1dXAErKylmwbLPPEYlIkCmZRKiTTIbFZhNXjYSEBG6+eqJXfu2dDZRXVPoYkYgE\nmZJJWHV1NTsKDnvlMcP6+RhNy8y4aBQ9umUCofkx+e/v8DkiEQkqJZOwfYdOejsZ5mRn0Mvn9bha\nIikpkRvzandj/MuS9TjnmrhDRCQ6op5MzGy2mW0zsx1m9lAj1zxpZjvNbJ2ZTYo4/6yZHTazDfWu\nf9jM9pvZ2vBrdlvjrNkIC2DE4N6YWVsf2SFmThtDelpooEDh0VN8sGWvzxGJSBBFNZmYWQLwFHAt\nMA64w8zG1LtmDjDcOTcSuA/4WcTbvw7f25CfOOcmh18L2hrrnv21yWTYoJ5tfVyHyUhPYdb087zy\nPG3tKyI+iHbNZAqw0zlX4JyrAJ4H5ta7Zi7wHIBzbhXQ1cz6hMvLgBONPLtdqw67933iHQ8b1Ks9\nHx1118+4gIRwTWrzrsI6tSwRkY4Q7WQyANgXUd4fPtfUNQcauKYh94ebxX5pZl3bEmR1dTUf7Y9I\nJgM7T80EoGdOF6ZPGu6VX9USKyLSwZL8DuAc/Q/wPeecM7PvAz8BvtjQhY888oh3nJeXR15e3qeu\nOXDklDesNic7g+5dM9s/4iibe+UElq/dBcCKtbu5+8Zpvm/qJSKdQ35+Pvn5+W16RrSTyQFgcER5\nYPhc/WsGNXNNHc65yHacXwB/aezayGTSmMhmoWEDO1cTV43hg3sxdng/tuw+SLVzvP7ORr5w83S/\nwxKRTqD+P7QfffTRVj8j2s1cq4ERZpZrZinA7cC8etfMA+4BMLNpwEnn3OGI9416/SNmFjk9/RZg\nU1uC3BPRXzK0E3W+13fTVbXDhN98b6uWWBGRDhPVZOKcqwLuBxYBm4HnnXNbzew+M/ty+Jr5wEdm\ntgt4Gvhqzf1m9r/ACmCUme01s78Ov/UjM9tgZuuAGcCDbYkzciTX8E7W+R7ponG53hIrpWXlvL1y\nm88RiUhQRL3PJDxsd3S9c0/XK9/fyL13NnL+nnaMj4LC4155yIAe7fXoDmdm3HjlBJ5+YSkAr7+z\nkeuuOJ/ERM1NFZHoCvy3zImiEs6UngUgLTW5U8x8b0relFF0yUgF4OiJYlZu+MjniEQkCAKfTPYd\nqp3GMqhvTqeZ+d6YlOQkZl9+vlfWXici0hECn0z2RjRxDe7X3cdI2s+1l471mra2f3SI3Xs1iVFE\nokvJ5GD8JZPuXTO5ZGLtJMbXl6p2IiLRFfhksu9QbTIZFCfJBOD6GbVNXcvW7uJkcYmP0YhIvAt0\nMnHOsfdgbZ9JvNRMAEbm9mFkbm8AqqqqWbR8i88RiUg8C3QyOXriNGfLQ3uYdMlIpVtWus8Rta8b\nZoz3jhcu20JlZZWP0YhIPAt0MqnfX9LZR3LVN23CUHKyMwA4WVzCinW7fY5IROJVoJPJvohkMqhv\n/DRx1UhKSuTay8Z55dfyN2onRhGJikAnk8Ijp7zjAX26+RhJ9FxzyViSkhIB2L3vKDsLjvgckYjE\no0Ank4NHa5NJ/97xmUy6ZqVz+YUjvPJrmsQoIlGgZBLWr1eb9teKaddfcYF3/N66PRw7edrHaEQk\nHgU2mZSUlntzLxITE+jdvXOvydWUoQN7MnZ4PyC0q+TCZRomLCLtK7DJpE6tpGdXEhLi+1dxXUTt\nZNGKLd7OkiIi7SG+v0GbEJQmrhpTLhjibeNbfKaMZR/s8jkiEYkngU0mhUdPesdBSCaJiQnMjhwm\n/I6GCYtI+wlsMqk7kiv+kwnA1dPPIzk8TLig8Bhbdh/0OSIRiReBTSaRc0yCUDMByMpMI2/KKK88\nX8OERaSdBDKZOOcoPFLbzBWvc0waEtkRv2rDRxw5XuxjNCISLwKZTIrPlFFSVg5Aakqyt35VEAzu\n150LRg0AwAEL3t3kb0AiEhcCmUwOHyvyjvv0yIq7BR6bc/2M2trJmyu2Una2wsdoRCQeBDKZHDle\nOwO8T49sHyPxx0XjcunbM/Rzl5SV887qHT5HJCKdXTCTSUTNpHePLB8j8YeZMfuy2p0Y5y/dpGHC\nItImgUwmRyNqJr1ygpdMAK6aNprUlGQA9h8+wfrt+32OSEQ6s0AmkyPHg10zAchMT2XmtNFe+XUN\nExaRNghkMomsmfTuHsxkAjDn8vOpGXqwdsteDkQMlxYRaY3AJRPnXJ3RXEGtmUBofs3ksble+Y2l\nGiYsIucmcMmk6HQZFZVVAGSkpZCZnupzRP66Pq92mPDiVds5U3rWx2hEpLMKXDKJ7C/pFeAmrhrj\nRw1gYJ8cAM6WV7B45XafIxKRziiAyUT9JZHMrM4kxvlLN1JdXe1jRCLSGQUvmdSZ/R68CYsNmXHx\nSK+578jxYtZs3utzRCLS2QQvmUQsbNgrjrfqbY3UlGRmXXKeV379nQ0+RiMinVHgksnRiGTSWzUT\nz+zLxpEQXqNs085CCgqP+RyRiHQmAUwmkbPfVTOp0at7FlPGD/XKmsQoIq0RuGRy7NQZ77hHt0wf\nI4k9N0R0xC9ds5Oi06U+RiMinUmgkknZ2QpKw/uYJCYmkJWZ5nNEsWXMsL4MHdgTgIrKKt58b6vP\nEYlIZxGoZBJZK+menRm4fUyaY2Z1aicL3t1MZXiCp4hIUwKVTI6fjEgmauJq0KWTRtA1Kx2A46fO\nsHLDRz5HJCKdQbCSSWTNpKuSSUOSkxO55tKxXlkd8SLSEoFKJsciaiY9lEwade2l40hMDP3V2PHx\nYXZ8fNjniEQk1gUqmZwoUjNXS+RkZ3DZ5BFeed4STWIUkaYFKpkcV82kxW66crx3vHLd7jrL9ouI\n1BeoZBI5miuna4aPkcS+IQN6Mn7UQAAc6jsRkaYFKpmoA751boyonbz13jbtdSIijWoymZjZbyKO\nvxD1aKLIOceJotoZ3Zr93rxJ5w1iUN/avU4WLd/ic0QiEquaq5lMiDh+IJqBRNup06XePh2Z6amk\nJCf5HFHsM7M6tZP5SzdpEqOINKi5ZOI6JIoOoAmL5+aKC0fVmcS4/MPdPkckIrGouWQy0MyeNLP/\nG3HsvVryAWY228y2mdkOM3uokWueNLOdZrbOzCZFnH/WzA6b2YZ61+eY2SIz225mC82sa3Nx1Fng\nUf0lLZacnMicy8/3yvOWbMC5uPk3hoi0k+aSyT8BHwBrIo4jX00yswTgKeBaYBxwh5mNqXfNHGC4\nc24kcB/ws4i3fx2+t75vAW8550YDi4FvNxdLZM1EI7laZ/Zl40hOSgTg4wOfsGlnoc8RiUisabLj\nwDn32zY+fwqw0zlXAGBmzwNzgW0R18wFngt/3ioz62pmfZxzh51zy8wst4HnzgVmhI9/C+QTSjCN\nOl5U4h13z1bNpDWyMtO4auoYFi7fDMC8Jeu5YNQAn6MSkVjSZDIxs3lNve+cu6mZ5w8A9kWU9xNK\nME1dcyB8rqk1PHo75w6HYzhkZr2biaPO3hzdstObu1zquSHvAhYt34wD1m7Zy96Dxxncr7vfYYlI\njGhuSNN0Ql/0fwBWAbG6ZnujjfiPPPIIAMvX7uJscm96DRxF1yw1c7VW/97duPiCIby/8WMAXl28\nnq/fdaW/QYlIu8jPzyc/P79Nz2gumfQFZgF3AHcCrwN/cM5tbuHzDwCDI8oDw+fqXzOomWvqO1zT\nFGZmfYEjjV1Yk0y++8QrbNtzCIBsbYp1Tv7q6kleMlm6Zie3z7mIXt2z/A1KRNosLy+PvLw8r/zo\no4+2+hlNdsA756qccwucc18ApgG7gHwzu7+Fz18NjDCzXDNLAW4H6jedzQPuATCzacDJmiasMOPT\nNaJ5wL3h4y8ArzYXSFFxbTNXzVBXaZ1RQ/owdng/AKqrq/lLvhaAFJGQZpdTMbNUM7sF+H/A14An\ngT+35OHOuSrgfmARsBl43jm31czuM7Mvh6+ZD3xkZruAp4GvRnz2/wIrgFFmttfM/jr81uPALDPb\nDswEHmsullOny7zjrl2UTM7VX13tjdzmzRVbtU+8iADNd8A/B5wPzAcedc5tau0HOOcWAKPrnXu6\nXrnBmo5z7s5Gzh8Hrm5pDJWVVd66UgZkZaa29FapZ9J5g8jt34OCwmOUV1Qy/91N3D7nYr/DEhGf\nNVcz+TwwktBSKu+ZWVH4VWxmnWZN8qIztbWSrC7pJCQEan3LdmVm3BJRO5n/zibKzlb4GJGIxILm\n+kwSnHNZEa/s8CvLOZfdUUG2VWRTjDrf2276xGH06RH6z3+m9Cxvrtjqc0Qi4rfmVg1OM7O/N7On\nzOzLZtYpV0es01+izvc2S0xMYO5VtWuAzluyXgtAigRcc+09vwUuAjYC1wE/jnpEURA5kitbne/t\n4sqpo+ssALl0zU6fIxIRPzWXTMY65z4f7jD/LHB5B8TU7k5GDgvuomau9pCSnMQNM2qXp//zWx9q\nAUiRAGsumXg9q865yijHEjXFZ9TMFQ3XXjaW9LQUAAqPnuK99Xt8jkhE/NLs5liRI7iA8Z1xNNfJ\n4tpFHjXHpP1kpqcy+9KxXvnFBR+odiISUM2N5kqsN4IrqXOO5qqtmajPpH3deOUEb9fKvQePs2rD\nRz5HJCJ+CMSEi1OntZRKtHTNSmf2ZeO88osL16p2IhJAgUgmdeaZqAO+3c2dOaHO5lmrNxX4HJGI\ndLRAJBOtyxVd3bIy6tVO1HciEjRxn0wqKqooLSsHICEhgS4ZWpcrGubOnOjVTvbsO8oHW/b6HJGI\ndKS4Tyan6i2lYhar+3t1bjnZGVwTMbLrhTfWqHYiEiBxn0wi55iovyS6bp45kaRw7WT3vqN8uHVf\nM3eISLyI+2RyuuSsd6wmrujq3jWTay45zyu/sEC1E5GgiPtkUrOPCYQm2Ul0zb1qIomJob9WOwuO\nqHYiEhBxn0xKSsu940zVTKKuZ04Xrp5WWzv539ffV+1EJADiPpmciUwm6Sk+RhIcn712sjey66P9\nn7BindbsEol3cZ9MTkc0c2UomXSI7l0zuX7GBV75+dffp6qq2seIRCTa4j6ZlET2maSpmauj3Dxz\nIhkRKwrnr97uc0QiEk1xn0wim7k0mqvjZGWmMXfmRK/8xzfWUF7RaXcxEJFmxH8yKVEzl19umHGB\nt0rzsZNnWLhsi88RiUi0xH8yUQe8b9JSk/nsNZO98stvrvWWthGR+BKAZKJ5Jn665pKx9MrJAkKr\nEcxbssHniEQkGoKVTNRn0uGSkxO5bc5FXnnekvWcKi5t4g4R6YwCkEzUzOW3GRePZGCfHADKzlbw\n4sIPfI5IRNpb3CeTsrMVABh4Q1WlYyUkJHDXjVO98sJlm9l36ISPEYlIe4v7ZFIjPS1Fy8/76OLz\nczl/ZH8Aqp3juVff8zkiEWlPgUkm6nz3l5lx782XUJPO127Zy7ptWgRSJF4EJplojon/hg7syZVT\nx3jl3/x5hZZZEYkTgUkm6nyPDXfeMIXUlGQA9h06wVvvbfU5IhFpDwFKJmrmigU52RncMmuSV37+\njTV1hm+LSOcUnGSiOSYx46Yrx9MzpwsARadLeXnRWp8jEpG2Ck4yUTNXzEhJTuLuG6d55dfe2cjB\no6d8jEhE2iowyUQd8LHl0snDGTWkDwBVVdX8+k8rtCOjSCcWmGSivUxii5nxxVsu9YYKf7ClgNWb\nCnyNSUTOXWCSifYyiT0jcntz9SW1+8X/6uXlnC2v8DEiETlXgUkmauaKTXfdMJWszDQAjp4o5uVF\nH/ockYici8AkE3XAx6aszDTuvql23a5XFq/jwJGTPkYkIuciQMlEzVyx6qqpY+p0xj/zwlJ1xot0\nMoFJJulaMThmmRn33Xo5CeGFODftLGTJqu0+RyUirRGYZJKSnOh3CNKEIQN6cuOV473yb155jxNF\nJT5GJCKtEZhkkpQYmB+107ptzkX06ZENhHbIfPbl5T5HJCItFZhv2KRE1UxiXWpKMl+57Qqv/N66\n3by/8WP/AhKRFgtMMklOUjLpDMaPHkjelNFe+RcvvquFIEU6gcAkk6SkwPyond69N08nu0s6AMdP\nneFXf1rhc0Qi0pxAfMMmJiZoy95OJCszjS9/7nKvnP/+djV3icS4qCcTM5ttZtvMbIeZPdTINU+a\n2U4zW2dmE5u718weNrP9ZrY2/JrdVAzqL+l8pk8cxmUXjvDKP3v+HU4Vl/oYkYg0JarJxMwSgKeA\na4FxwB1mNqbeNXOA4c65kcB9wM9beO9PnHOTw68FTcWRrCauTulLn72cnOwMILTviSYzisSuaH/L\nTgF2OucKnHMVwPPA3HrXzAWeA3DOrQK6mlmfFtzb4nYr1Uw6py4ZqXz1jjyvvHLDRyxds9O/gESk\nUdFOJgOAfRHl/eFzLbmmuXvvDzeL/dLMujYVhEZydV6Txw7m6um1Kws/8+K72khLJAYl+R1AA1pS\n4/gf4HvOOWdm3wd+AnyxoQu3rHyNrMxUHnlkJ3l5eeTl5bVjqNIR7r15Opt2HuDQJ0WUna3gv557\nmx88MJck/SNBpF3k5+eTn5/fpmdEO5kcAAZHlAeGz9W/ZlAD16Q0dq9z7mjE+V8Af2ksgLHTbmBQ\nv+488q1bWx28xIb0tBQevOdqvvPEK1RVVbNr7xGen7+az980rfmbRaRZ9f+h/eijj7b6GdFu5loN\njDCzXDNLAW4H5tW7Zh5wD4CZTQNOOucON3WvmfWNuP8WYFNTQWgplc5vRG5v7rx+ild+5e11rN++\n38eIRCRSVL9lnXNVwP3AImAz8LxzbquZ3WdmXw5fMx/4yMx2AU8DX23q3vCjf2RmG8xsHTADeLCp\nONRnEh/mXjWB8aMGAuCAJ373NsdPnfE3KBEBwOJ5qKWZuVu+8TPGDu/Hv3+j/iAy6YxOFJXwzcdf\npOh0aM7JmGF9efRrN6r/RKQdmRnOuVbN9A5E+49qJvEjJzuDB++Z6Y3S2LbnEL9/7X1fYxKRgCQT\nzTOJL+NHD+SOG2r7T+YtWc+Kdbt9jEhEApFMNAM+/txy9SQuGpfrlZ/6fT57Dx73MSKRYAvEt2yi\nmrnijpnx9c9f5W2mdba8gsd+sYDiM2U+RyYSTIFIJuoziU9dMlJ56G+vJTUlGYDDx4r4j18torKy\nyufIRIInIMkkED9mIOX278EDd1/llTfvKuSXLy/TgpAiHSwQ37KqmcS3qeOHckfEhMY3V2zltfyN\nPkYkEjyBSCYazRX/PjNrEpdOrt3/5LevrGD5hxrhJdJRApFMVDOJf2bG/XfmMXpoaKWdmhnym3cV\n+huYSEAEIpkkam2uQEhJTuLbX5pN/16hHQmqqqp57BcLNGRYpAME4ltWNZPgyMpM41+/egPdskI7\nNJaUlfO9/3mNQ58U+RyZSHwLRDJJ0miuQOndPYt/+cp13pDhE0UlPPLUX/jkxGmfIxOJX4H4llXN\nJHiGDuzJt7802/tvf/REMY/+9184WVzic2Qi8SkQyUT7mQTTBaMG8E9/c43XZ1Z49BSP/vdrnCou\n9TkykfgTiG9Z1UyC68JxuTx4z9XeKsN7Dx7n4afmcaJINRSR9hSIZKJ5JsE2feIwvv75q7yEsu/Q\nCf7tyVc5dlJ9KCLtJRjJRB3wgTfj4lH8/T1Xk2ChlFJ49BT/+uQ8jhwv9jkykfgQiG9Z7cInAJdd\nOIJv3jvL60M5fKyI7/z0z3x84BOfIxPp/AKRTNRnIjWmTxzGP3/xWi+hnCgq4V+enMfGHQd8jkyk\ncwtEMtFoLol00bhc/u3vric9LQWA0rJy/v3nr7Psg10+RybSeQXiW1Y1E6nv/JED+MEDc8nJDs2U\nr6qq5qfPvcUf31ij5etFzoGSiQRWbv8e/PDBv2JA727euRcWrOE/f/0mZWcrfIxMpPMJRDLRQo/S\nmN7ds/jhg3/F+FEDvXMr1+/hu0+8qpFeIq0QiG9Z1UykKV0yUvmXr1zHdVec7537+MAn/MPjL7J6\n08f+BSbSiQQimagDXpqTmJjAFz9zGffdegUJCaG/LyVl5Tz2iwX8bt5KqqqqfY5QJLYF4ls2OVk1\nE2mZay4dyw8emEuPbpneuVfeXsd3n3iFg0dP+RiZSGwLRDLRcirSGqOG9OHH//w5Jp03yDu3s+AI\n//Cjl3g4LR1rAAAMTklEQVTrva0a7SXSgEAkk2QtpyKtlJWZxnfvu467bpjqNXudLa/gZ8+/w2O/\nWKB1vUTqCcS3rGomci7MjFtmTeKxesOH12wu4IH/8wILl21WLUUkLO6TiaGhwdI2wwf34j/+6TPM\nubx2tFdpWTnPvPgu//rkPAoKj/kYnUhsiPtvWS3yKO0hNSWZv/3sZXzv6zfRv1dX7/zWPQf5xx+9\nxLMvL+N0yVkfIxTxV9wnE80xkfY0bkR/fvzQ5/jMrMleX0q1c8xfuon7v/8H3nh3E5WVVT5HKdLx\n4j6ZqGYi7S0lOYk7b5jCj//5s5w/sr93vvhMGb98aRn3f/958t/fTnW15qZIcFg8dyCamfvSvz3H\nM4/e7XcoEqecc6xc/xG/eWUFn5yoO8JrYJ8c7rj+YqaOH4qZNfIEkdhjZjjnWvWXNilawcQKjeSS\naDIzpk8cxoXjBvPGu5v505trvb6T/YdP8B+/WsTgft25eeZELp00XDVliVtxXzP5xg+e54nv3OZ3\nKBIQJaXlzMtfz7zFGzhbXnfl4R7dMrl+xnhmTT+PjPQUnyIUad651EziPpk8+NgL/OShz/kdigRM\n0elS/vTmhyxcvoXyiso676WnpXDV1NFcPf08Bvfr7lOEIo1TMqnHzNw//+fLPP4Pt/gdigRU8Zky\nFizbzPylmyg6Xfqp90cN6cM1l4zlkknDSE1J9iFCkU9TMqnHzNx3/uvP/OCBm/0ORQKuvKKSd1bv\nYN7i9RQ2sGBkeloK0ycM47ILR3D+iP6aaCu+UjKpx8zcw0/N45Gv3eh3KCJAaPTXhh0HWLR8C+9v\n/LjB4cPZXdKZPmEYl04eztjh/TQSTDqcRnM1QHuZSCwxMyaMHsiE0QM5VVzKkve389Z7W+ssb190\nupSFyzezcPlmsrukM3HMQC4cl8uk8waRmZ7qY/QijYv7msljv3iDh/52tt+hiDTKOceOjw+z/MPd\nrPhwNyeKShq8LsGMMcP6MmHMIM4f0Z8Rg3tpqLFEhZq56jEz9x+/WsQ//vUsv0MRaZHq6mq27jnE\n8rW7WblhD6eKP91pXyMlOYnzhvVl3Mj+jBnal+GDepGWqk58aTslk3rMzP3Xc2/xwN0z/Q5FpNWc\nc+zee5Q1WwpYu3kvu/cdbfL6BDMG9evOqCG9GZnbm+GDejGgd452GpVWUzKpx8zcU79fwtfuzPM7\nFJE2O1FUwvpt+9i0q5DNOws5cry42XsSzOjfuxu5A3owuF93cvt3J7d/D3rldFHHvjRKyaQeM3NP\n/3EpX771cr9DEWl3R44Xs3lnIZt3F7Lz4yMcOHyClv7fnJSUSN8e2fTr1dV79e2ZTd9eXemenaG+\nmICLydFcZjYb+C9CKxQ/65x7vIFrngTmAGeAe51z65q618xygD8CucDHwK3OuU8P3kdL0Ev86t09\ni95TR3Pl1NFAaCmX3fuOsqPgMLsKjlBQeJzDx4oavLeysor9h0+w//CJT71nQNesDHp0y6RHt0x6\n5nShe9dMenbrQvdumXTLziA7M40uGamq3YgnqjUTM0sAdgAzgUJgNXC7c25bxDVzgPudc9eb2VTg\nCefctKbuNbPHgWPOuR+Z2UNAjnPuWw18vvvdq+/x+ZumRe1n7Czy8/PJy8vzO4yYEKTfRWlZOfsO\nnaCg8BgFhccpKDzGvkMnKD5TBsDR/TvoNXDUOT07wYysLmlkd0mna5c0sjJDf3bJTCMzPYXM9BTS\n01LITE8lIy2ZjPRUMtJC51OSY29WQpD+XjQnFmsmU4CdzrkCADN7HpgLbIu4Zi7wHIBzbpWZdTWz\nPsDQJu6dC8wI3/9bIB/4VDIBSFTNBND/KJGC9LtIT0th1JA+jBrSp8750yVnOXT0FD/8wfeZPeci\nDh49xaFPijhyrJhTxSUtai6rdo5TxaWcKi5lXyvjSkxMID01mbTUZFKTk0hJSSI1JYnU5NCfKRHH\nNe+nJCeRlJhAUmIiycmhPxMTE0LnkhJJTkokKTHB+zMxMZGkpFA5wYyEBCMxIYGEhNBxzbkESyAx\n0ViyZElg/l5EQ7STyQCo8/dsP6EE09w1A5q5t49z7jCAc+6QmfVuLAA1c4l8WpeMVEbk9mZw/+7c\nOvuiOu9VVlZxvKiE4yfP8MnJ0xw7eSbi+DRFp8soOlNGaVn5OX9+VVU1p0vOxtRWx1tWfsDmkz8n\nITEhnGgS6iadegkIwCz0r/ia1j4zo+af86Hz5l1HE+/V3PepZ0VeQ/33Pv0zRDY71n5a/Wsaub6N\nLZaxV9ekkd9A0xr9h5RmwIu0TlJSYqg/pntWk9dVVFRRdKbUSy5FxaWcOl1KcclZSsvKOVNaHv7z\nLCVlFZSUnuVMaTklZeVUVcXmLpSOUKILbbys7ZdbxTkXtRcwDVgQUf4W8FC9a34O3BZR3gb0aepe\nYCuh2glAX2BrI5/v9NJLL730av2rtd/30a6ZrAZGmFkucBC4Hbij3jXzgK8BfzSzacBJ59xhM/uk\niXvnAfcCjwNfAF5t6MNb24EkIiLnJqrJxDlXZWb3A4uoHd671czuC73tnnHOzTez68xsF6GhwX/d\n1L3hRz8OvGBmfwMUALdG8+cQEZGmxfWkRRER6Rhx2TttZrPNbJuZ7QjPQwkkMxtoZovNbLOZbTSz\nb/gdk9/MLMHM1prZPL9j8VN4CP6LZrY1/Pdjqt8x+cXMHjSzTWa2wcx+b2YpfsfUkczsWTM7bGYb\nIs7lmNkiM9tuZgvNrGtzz4m7ZBKe7PgUcC0wDrjDzMb4G5VvKoFvOufGAdOBrwX4d1HjAWCL30HE\ngCeA+c6584AJhAa1BI6Z9Qe+Dkx2zo0n1PR/u79RdbhfE/q+jPQt4C3n3GhgMfDt5h4Sd8mEiImS\nzrkKoGayY+A45w7VLE3jnDtN6AtjgL9R+cfMBgLXAb/0OxY/mVk2cLlz7tcAzrlK51zD664EQyKQ\naWZJQAahFTcCwzm3DKi/rs5cQhPCCf/Z7N7n8ZhMGpsEGWhmNgSYCKzyNxJf/RT4J0JDH4NsKPCJ\nmf063OT3jJml+x2UH5xzhcCPgb3AAUKjSd/yN6qY0DtyYjjQ6MTwGvGYTKQeM+sCvAQ8EK6hBI6Z\nXQ8cDtfUjHObHBsvkoDJwH875yYDJTSyHFG8M7NuhP4Vngv0B7qY2Z3+RhWTmv0HWDwmkwPA4Ijy\nwPC5QApX3V8Cfueca3A+TkBcCtxkZnuAPwBXmtlzPsfkl/3APufcmnD5JULJJYiuBvY4544756qA\nPwGX+BxTLDgcXiMRM+sLHGnuhnhMJt5EyfCojNsJTXIMql8BW5xzT/gdiJ+cc99xzg12zg0j9Hdi\nsXPuHr/j8kO4+WKfmdUsFzyT4A5K2AtMM7M0Cy1UNZNgDkaoX1uvmRgOTUwMjxSLa3O1STOTHQPF\nzC4F7gI2mtmHhKqq33HOLfA3MokB3wB+b2bJwB7Ck4WDxjn3vpm9BHwIVIT/fMbfqDqWmf0vkAf0\nMLO9wMPAY8CLrZkYrkmLIiLSZvHYzCUiIh1MyURERNpMyURERNpMyURERNpMyURERNpMyURERNpM\nyUSkhcysuIFz95nZ58PHo83sQzP7wMyGNvGcb9crL2v/aEU6luaZiLSQmRU557KbeP8hINE598Nm\nnlPsnMtq9wBFfBR3M+BFOpKZPQycJrQcyd8DlWY20zk308zuIjTTPJnQas1fA34ApJvZWmCzc+7u\nmuRiZjOAR4GTwPnAi8BGQnuwpAE3O+c+MrOewM+BQeEwHnTOreion1mkIWrmEmk755x7g9AX/E/D\niWQMcBtwSXhl3mrgTufct4ES59xk59zdNfdHPGs88GVgLHA3MNI5NxV4ltAmThDa2Oon4fOfJeD7\ns0hsUM1EJDpmElqJd3V4AcE04FD4vaaWv1/tnDsCYGa7Ca0xB6EaSl74+GrgvPBzIbRseoZzrqQd\n4xdpFSUTkegw4LfOue+28r6zEcfVEeVqav9/NWBqeCdRkZigZi6RlmvNhlpvA581s14AZpZjZjV9\nHOXhfWbO5bkQqq084N1sNqGV94u0OyUTkZZLN7O9ZrYv/Off08gOdOFtD/4FWGRm6wklgH7ht58B\nNpjZ72oub+TzGjv/AHCRma03s03Afefyw4i0Jw0NFhGRNlPNRERE2kzJRERE2kzJRERE2kzJRERE\n2kzJRERE2kzJRERE2kzJRERE2kzJRERE2uz/A00f4OrL9mK+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f48daeb8510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf = MakeWeibullPmf(lam, k, high=10)\n",
    "thinkplot.Pdf(pmf)\n",
    "thinkplot.Config(xlabel='Lifetime',\n",
    "                 ylabel='PMF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use np.random.weibull to generate random values from a Weibull distribution with given parameters.\n",
    "\n",
    "To check that it is correct, I generate a large sample and compare its CDF to the analytic CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlwnPWd5/H3tw+1Wmp1S63Lso1tbrLUkoQAQ47NKsns\n4mRrIDO1M4GkspBkUuwWJOxRGzKZmopmajaZpCZhwmASHBwCwxkwMeayjQ8ZbGzj+z7kW75tLINv\nS+rf/tHNY0m2rKulp4/Pq8qV53n66ae/dOyPfvo+z/N7zDmHiIgUloDfBYiISPYp3EVECpDCXUSk\nACncRUQKkMJdRKQAKdxFRApQn+FuZlPM7KCZrbnEPg+bWYuZrTKzT2S3RBERGaj+jNyfAG7r7UUz\n+zJwpXPuauBe4DdZqk1ERAapz3B3zi0A2i6xyx3AU5l9lwAJM6vPTnkiIjIY2ei5jwFau6zvzWwT\nERGf6ISqiEgBCmXhGHuBy7qsj81su4CZaSIbEZFBcM7ZQPbvb7hb5s/FTAfuA14ws1uBY865g5co\ncCD1Faympiaampr8LiMn5NJ30d7eyb7DxzjSdoKW3YdwKccHJ06zacdBWvcfBeAYUdoC5aR6/Scx\neAcWvcCoT39twO8rLQlhGKfPtVNRFsHMcM5xrr2TivIIp8+0EytPbzfALF17e0cn0UgYzAgYgNGZ\nSlESDpJeM+9f/kfLH/1Xp5wjFLzwl3/r+b1c5Gu6YJNduNPaN6Zww1e+08sb7eKbL8IC518NBwME\nAvnXsJjy3ZsH/J4+w93MngUagWoz2w38GCgBnHNusnPuDTP7ipltBU4C3xpwFSIjrKOjk627D9O8\ndDPL1+/mxKmzdKZSvN8Z6SW4gxCozcpnx8oiBMw4195BOBykPBrBpRyRSJgzFVEuG1VFOBQkGAjg\ncIRDQQJmYEZZJMiffWI0f/apMV5AF6qmI7Nouv8zfpeRE6Z8d+Dv6TPcnXNf78c+9w/8o0VGxsnT\nZ9m6+zB7D7axdO0uFm5r44iLdgnwsswfBnQWKhoJEwoFCAWDtLd3UB6NpI8WLSEUChIKBDBLj6q/\n+qkx3H7T2D6P2XRqPk0PfnFg/4EiF5GNnrsMQmNjo98l5IxsfBenz5xj78Fj7D10jN37j7J+6z5a\ndh3qpYVSdsljBcxIOUdVvAwzKAmHiEbCBIMBSsIhSktCRCMhbv/kaG67YdSQa+9Kfy/O03cxNDaS\nPXAzc+q5y1CdOdvOolXbWbZ+F8FggIUrtgKD74XHyiIkYlESsVKipSWUlgSHJbhFBitzDmVYTqiK\n+OrQ0eO8vayF515/r9v2dKDX9DvQqyqixMpLqY5HufOzE5h4Q8NwlCviO4W75KSTp8+yckMrG7fv\nZ9a7G0mlUt5rfY3QA2ZUV5YTCgaIx6JUVZTy5zeN1UhciorCXXLGlp0HmT5vDfsPf8DOvUe87b2F\neTQSpixaQnWinPJohGAwQCQcUEtFBIW7+Oz4yTP82/TFzFm8qdv23gJ9VE2c6spyopESb5sCXeRC\nCnfxxa597/PG2+uYs2gjXU+xdw31aCRMaSRMdWU58fLSbjefKNBFLk3hLiPGOceb76xjxjvr2Xvo\nmLe9a6AHA8ao6jhV8TKipRqdiwyWwl2GXWdnikWrtvPQU7O7bf8o1APBIKNrElRXlhMOBb3XFegi\ng6dwl2G1Y88RHnm22TtB2nWUXloSYmx1nLpkjK6zgyjURYZO4S7DoqOjk6lvreQPM5YBF54gvXp8\nHYlY1NtfgS6SXQp3ybr5S7fw3OtLOdx2vFuo1yUrqK+uIFIS9vZVqIsMD4W7ZM2m7Qf43csL2dZ6\nGEiP1t8PxCgtCTF+dDUV5aWAAl1kJCjcZchSqRQvzFjOSzOXA91bMPXJCsaOqsLMFOoiI0jhLkNy\n9IOT/PS3M9jeerhbqJeXljB2VJU3Wv/LW3T7v8hIUrjLoL23dieTnp3HiVNnvRZMKBhg/Kgqqitj\ngFowIn5RuMuAHT56nEefm8+aLXu6zcpYWxVjXENSLRiRHKBwlwGZu3gTk55r7hbqwYBx9dgaEhXp\nh2CoBSPiP4W79Mvxk2eY/OI7vLtym9eCAaisiDKuIUlJOKTRukgOUbhLn2YuWM/kF9+54MEY4xqq\nqEvGFeoiOUjhLr3q6Ojkl0/OZuaaA91CPRErZcKYasKhkFowIjlK4S4XdeZsOw89OZtl63d1C/a6\nZAXjGqqIhPWcUZFcpnCXC7x/7AR/P+k19h46xjGipDBKQkGumVBPIhZRqIvkAYW7dLNs/S5+OvnN\nbv31qngZV4ytobQkyKS7b/S7RBHpB4W7eF6Zu5qHX1nVrQ0zqibO2PpKwLj9k6P9LVBE+k3hLjjn\nePrVJUyb0z3YLxtVRX21roYRyUcKd2HqWyuZNmeV11+PRsJcNa6WeLn66yL5SuFe5KbPW82vX1/r\njdjLS0u4ZkI9ZaUh9ddF8lig712kUL0+fy1PTlvUbSbHay+vJxgMqL8ukuc0ci9SC5Zv5Zcvr/BG\n7NFImKvH1xGNhNSKESkACvcitG33YR59fj5tgQQpjEg4yLUT6imPhtWKESkQCvcis2pTKz98fD6H\nOhOZGR0DfOyKBsqjYbViRAqIwr2IbN11iB9MbuaQS0/NGzDjqvG1GrGLFCCFe5FYvn4XP5sykyMu\nCUDA4KpxtdQkyjRiFylACvcisGvf+/xgcjNtgSQpjIDBdVc0cPfnr9CJU5ECpXAvcPsOHaNp0mve\n5Y4BM669vJ6qilIFu0gB69d17mY20cw2mdkWM3vwIq/HzWy6ma0ys7Vmdk/WK5UBO3n6LN97+C1W\nnyr3phS4ekIdyXhUrRiRAtfnyN3MAsAjwJeAfcBSM3vFObepy273Aeudc7ebWQ2w2cyeds51DEvV\n0ifnHD+fMpMdJ80L9isvq6EmUaaTpyJFoD8j91uAFufcLudcO/A8cEePfRxQkVmuAN5XsPvHOcdT\nryxmQUubF+zjG5KMqq7QiF2kSPSn5z4GaO2yvod04Hf1CDDdzPYBMeBr2SlPBmPqWyuZPm81bYEa\nAJKJcsbWJzRiFyki2Tqhehuw0jn3RTO7EnjLzG5wzp3ouWNTU5O33NjYSGNjY5ZKEICVG1u7TQRW\nVlrC5WOqNWIXySPNzc00NzcP6RjmnLv0Dma3Ak3OuYmZ9R8Czjn3sy77vAb81Dm3MLM+B3jQObes\nx7FcX58ng3f0g5M88JMXWH8u7j0a7/qrRmuGR5E8Z2Y452wg7+lPz30pcJWZjTezEuBOYHqPfXYB\nf5opoh64Btg+kEJkaJxzPPDIbC/YQ8EA110xirLSkEbtIkWoz7aMc67TzO4HZpH+YTDFObfRzO5N\nv+wmA/8I/N7M1mTe9gPn3NFhq1ouMHfJJtYfafdOoF4+poaKshKN2EWKVJ9tmax+mNoyw+Lw0eN8\n4x+neXPG1FbFuGZ8rabuFSkQg2nL6A7VPJdKpfjlk7M54qIAhIIBrhxbrRG7SJFTuOe5v3n8bd7a\n3dmtHfPnN431uSoR8Zses5fHlqzZwVsbj3jBXp+soC5ZrlaMiGjknq/+8O4Ofj51hRfssWiEKy/T\n9ewikqZwz0MdHZ08/OpaOlLp9WDA+NgV9Tx6z6f8LUxEcobaMnnoH55ZzIkz7d76dRPq+K+3XOZj\nRSKSazRyzzNL1+3ktVX7IdOOGVef4Jnvfc7fokQk52jknkc+PHGax154u1uf/b9P/Hc+VyUiuUjh\nnkf+7smFrDiRvlEpGAhw5bhavvIJnUAVkQsp3PPEtt2HeWfrsfPzs49OEouGfa5KRHKVwj0PzFiz\nn+9MXuIFeyJWSkONHrwhIr3TCdU88Ls5Wzh28qy3Pr6hStMLiMglaeSe4w4c+ZBte9731huqY3zj\nc1f4WJGI5AON3HNYR0cnP/jtPDpT6Zk0g4EAU//PFyiLlvhcmYjkOo3cc9gzr73HukPnb1b62OV1\nCnYR6ReN3HPU7+ds5uH5e72TqKOq43zz81f6XJWI5AuN3HPQ2XPtTJ610Qv2stISrhxbrdkeRaTf\nFO45JpVK8dCTczh1rtPbdvW4Gm6/UZc9ikj/qS2TY37x8gr+sOEEH80dM7o2wZP/49P+FiUieUcj\n9xzinGPasvN99upEOZePTvpclYjkI4V7Dlm4chsnz3YAEDDjqsvUjhGRwVFbJke0t3fyzy8t56P/\nS+qSFTz2nZv8LUpE8pZG7jlievNqdp8OAuknK41vqPK5IhHJZwr3HLD/8Ae8OGO512tvqE3wFzeP\n9bkqEclnCvcc8Le/e4ctqfRIvbQkRH11XNe0i8iQKNx91rLrIKsOnPFG7RPGVFNaEvS5KhHJdwp3\nn73w5jIv2JOJcqoTZZqnXUSGTFfL+OihP67g5c2nwZs/pkLztItIVmjk7pPDR4/z7IKd3qi9Kl5G\nVTzqc1UiUigU7j75/bRFtGfmaS8JB7lmXI3aMSKSNWrL+ODp+S08t/ZDPmrHTBhdzWPfudnfokSk\noGjk7oMn5rR47Zh4eSm1VeU+VyQihUbhPsJ27j3C+8dPe+sTGqrUjhGRrFNbZoT95Ln3vOVErJSn\n7vuMj9WISKHSyH0Erd68h2V7T3nrExo0na+IDI9+hbuZTTSzTWa2xcwe7GWfRjNbaWbrzGxedsvM\nf845npy2yOu1J2JRvv65y32uSkQKVZ9tGTMLAI8AXwL2AUvN7BXn3KYu+ySAScB/ds7tNbOa4So4\nXz30x5U0H0gHuxmMH53U/DEiMmz6M3K/BWhxzu1yzrUDzwN39Njn68BU59xeAOfckeyWmd/OnG3n\nhUW7vFF7fXWcirISn6sSkULWn3AfA7R2Wd+T2dbVNUDSzOaZ2VIz+2a2CiwE0+au4mxHCsjM1T6q\nUlfIiMiwytbVMiHgRuCLQDmwyMwWOee2Zun4eWva0lZ+Mfv8z8aG2krdsCQiw64/4b4XGNdlfWxm\nW1d7gCPOuTPAGTN7G/g4cEG4NzU1ecuNjY00NjYOrOI8M2X2pm7TDIyrT/hckYjkuubmZpqbm4d0\nDHPOXXoHsyCwmfQJ1f3Ae8BdzrmNXfa5DvhXYCIQAZYAX3PObehxLNfX5xWSY8dP8cUfzyCV+W++\ndnwt3/7iNTqRKiIDYmY452wg7+lz5O6c6zSz+4FZpHv0U5xzG83s3vTLbrJzbpOZzQTWAJ3A5J7B\nXoz+4elFXrCXloR49nufIxjUrQUiMvz6HLln9cOKaOS+cdt+vvHIQu8Kmesm1PL8A5/3uSoRyUeD\nGblrGDlMnnzl/A1LsWiEe75wtc8ViUgxUbgPgy07D9Ky65C3fllDFRNvaPCxIhEpNgr3YfDzF5ex\nI5C+SbeyIkp5NOJzRSJSbBTuWbZjzxFW7DvttWRG1yaIhPU1i8jIUupk2dS3VnabHKwqHtXdqCIy\n4jSfexY9Pb+FZ9d8wEePz2uojTPp7hv9LUpEipJG7ln0VPO2bqP26kSZzxWJSLFSuGfJ/sMfcOjY\nSW99bH1C7RgR8Y3aMlny1rvnb8itKIvw9H2fwWxA9xyIiGSNRu5Z8OGJ0zy3YLu3XpusULCLiK8U\n7lkwbc4qDnakr2WPhIOMqo75XJGIFDuF+xCdOdvOnMWbzl/XXlfJHTf2fJaJiMjIUrgP0U+ef4+1\nZyqA9HztyUS5pvQVEd8p3Ifg5OmzvLnmoDdqH1WToLQk6HNVIiIK9yGZtXCD95SlcCjImNq4Ln8U\nkZygSyEH6Vx7B9PnrQHSJ0/rq+P8+luf8rcoEZEMjdwHadJra1h9qhyAUDBAXVJXyIhI7lC4D9K0\nZfu8XntdsoJoRL8EiUjuULgPwooNu/ng1FlvfbR67SKSYzTcHIRX5q7ylmsqy3n8u7f4WI2IyIU0\nch+g1Zv3sKClDUhP7Du6rtLfgkRELkLhPgCpVIqnX11CWyB9IrUqUU5FWYnPVYmIXEjhPgCL1+xg\ne+thUhgGjK2vVK9dRHKSwr2fOjo6eeGNZd56fXWcknBIUw2ISE7SCdV+emf5VtYdPENboAYzqK+p\n8LskEZFeaeTeT6/MW01boJwURn11nHAoRCSsr09EcpPSqR/Wb91H6/6j5ycIq44TCQfUbxeRnKW2\nTD88/8ZSbzmZKCMUCjLp7ht9rEhE5NIU7n3YufcI7247lu61Aw21uq5dRHKf2jJ9mLFgvddrr6yI\nEo2E1WsXkZynlLqEk6fPMn9pi9drr6mqUK9dRPKC2jKXsGD5Vg61hyEApSUhEhWl6rWLSF7QyP0S\nFqzY6k01kB616xF6IpIfFO69ONJ2gg3b9nstmWSiTO0YEckbCvdeLFy5zVuuKI9oqgERySvquV9E\nR0cnv3trI62BGgCq4mU+VyQiMjD9Grmb2UQz22RmW8zswUvsd7OZtZvZX2SvxJG3YmMrrWdCpDCC\ngQA1lTFd/igieaXPxDKzAPAIcBtwPXCXmV3Xy37/BMzMdpEj7dV5q71ee20yRjQSUr9dRPJKf9oy\ntwAtzrldAGb2PHAHsKnHft8DXgJuzmqFI2zn3iNs2LYfArUA1CYrdPmjiOSd/vQaxgCtXdb3ZLZ5\nzGw08FXn3K8hM+TNU795cz07Mr32RCxKJKzTEiKSf7LVSP4XoGsvPi8DvqOjk7mb3j/fkqlSr11E\n8lN/hqV7gXFd1sdmtnV1E/C8mRlQA3zZzNqdc9N7HqypqclbbmxspLGxcYAlD593lm/lXKcDIBwM\nUFdVrl67iIy45uZmmpubh3QMc85degezILAZ+BKwH3gPuMs5t7GX/Z8AXnXOvXyR11xfn+env/3V\nNF7f2QlAQ22cN3/0n3yuSEQEzAzn3IA6In2O3J1znWZ2PzCLdBtninNuo5ndm37ZTe75loEUkCt2\n7TvK4u0fQCAGQG2VHqMnIvmrX2cLnXMzgGt7bHusl32/nYW6RtyC5S3ePDKVFVEqykp8rkhEZPB0\ntpD0idRZ727wTqRWV6rXLiL5TeFO+kTqiVNnAQgGAlRWlGkeGRHJawp3YO6S8/dj1SVjpC/6ERHJ\nX0V/h86eg23eM1IhfUeqiEi+K/qR+9RZK7xnpFaUpaf21Y1LIpLvijrFzpxt591V270TqQ21CT0j\nVUQKQlG3ZVZtaqWjoxMCECkJEY/pGakiUhiKeuT+1LwWb5KwZLyMPJ0SR0TkAkUb7kfaTrBk93Gv\nJVMZL1OvXUQKRtGm2avz1tCZmaqhPFpCsiKqXruIFIyi7Ll/eOI0MxeuB5IA1FfHmXSPeu0iUjiK\ncuQ+f2kLhzvSc8dESkJ6ALaIFJyiC3fnHLMXbfQmCRtVE6e0JOhzVSIi2VV04b5my172HGwjhREw\nIxnXJGEiUniKLtzffHudt5xMlBEMBjRJmIgUnKI6ofr+sRPMXneQo5lr2+ur4z5XJCIyPIpq5D53\nyWaOZuaRiZVFiJaW6Np2ESlIRZNs59o7eK15TbcHcmgeGREpVEXTlml+b0v6gRwBKAkFqamMaR4Z\nESlYRTNyf2LuZm8emdpkhR7IISIFrSjCfcvOg2w+2kkKwwxqqmLqtYtIQSuKhJu7ZJPXa6+KlxOL\nhtVrF5GCVvA995Onz9L83hY+mkemprJcvXYRKXgFP3J/9PW1bElVARCNhInHSn2uSERk+BV8uL+6\ncn+Pyx81j4yIFL6CDvfd+4/y4alzAJjB6Nq4eu0iUhQKuuf+6rw13nIiFmXyd272sRoRkZFTsCP3\n4yfP8M7yFm+9LlnhYzUiIiOrYMN97pLN3R7IoROpIlJMCjLcnXPM6fJAjvpkhU6kikhRKchwX7mx\nlb2HjmUeyAHVlTGdSBWRolKQ4T7jnfXeclWiXA/kEJGiU3BXy7R9eIq5Gw55D+RoqEn4XJGIyMgr\nuJH7uyu3eQ/kqCiLUBoJa5IwESk6BZd685duOT9JWKJMD+QQkaJUUG2ZQ0ePs631MARqMYNkQpOE\niUhx6tfI3cwmmtkmM9tiZg9e5PWvm9nqzJ8FZvbvs19q3x57c733QI54eSmhoC5/FJHi1Ge4m1kA\neAS4DbgeuMvMruux23bg8865jwP/CPw224X2xTnHm2sOdpm3vUy9dhEpWv1Jv1uAFufcLudcO/A8\ncEfXHZxzi51zH2RWFwNjsltm37buPsSpcx0ABMwYVV2hXruIFK3+9NzHAK1d1veQDvze/DXw5lCK\nGozX56/zlqviZfzm2zeNdAkiIjkjqydUzewLwLeAz/W2T1NTk7fc2NhIY2PjkD/35OmzvL5qL5Ce\nbqCmKjbkY4qI+KW5uZnm5uYhHcOcc5fewexWoMk5NzGz/kPAOed+1mO/G4CpwETn3LZejuX6+rzB\neGdZCw88s5oURkkoyM3Xj2PSPbpKRkQKg5nhnLOBvKc/PfelwFVmNt7MSoA7gek9Pngc6WD/Zm/B\nPpxmLtxw/mlLVTFuv1G9dhEpbn22ZZxznWZ2PzCL9A+DKc65jWZ2b/plNxn4O9JPoH7UzAxod85d\nqi+fNVt2HmTj9v0QqAXSLRnNIyMixa5fPXfn3Azg2h7bHuuy/F3gu9ktrX/+9dXV3rXtyUQZkXBB\n3ZclIjIoeX0h+Okz51i047jXkqmvjuvadhER8jzcm5duoSNzfjZaEiZZEdW17SIi5PncMvOWbPaW\na5MxXSEjIpKRtyP3p+e3MHtvethuBsnKcp8rEhHJHXkb7s8t2On12isryigvDftckYhI7sjLcO/o\n6OTA0RPeeoPmkRER6SYve+6L1+ygM5UCIBwM8vT9nyEQyMufUyIiwyIvE/G3Mzd4yzXJcgW7iEgP\neZeKLbsOsvbQOQAMGFOrB2CLiPSUd+E+Y8GGbs9I/cs/GedzRSIiuSevwv3suXYWrtjqrdclKzSP\njIjIReRVuP/z1BVsSVUBEAmHiJVFfK5IRCQ35U24p1IpXlu132vJ1FTFiIT1AGwRkYvJm3BfsmYn\np891AhAwGFMb17XtIiK9yIvr3GeuOcBPp56//LEuGWfyX9/sY0UiIrktL0buT83fyrETZ731MXVx\nH6sREcl9OR/uzjl27Gvz1uuqyvmrW8f7WJGISO7L+bbMyo2tnDqTuWnJ4Lnvf47qypjPVYmI5Lac\nH7n/8o8rvOWaypiCXUSkH3I63DfvOMCGIx3e+viGKh+rERHJHzkd7i/OXO5d157UVAMiIv2Ws+G+\nrmUvKze2euujahKaakBEpJ9yNtx//uIydgRqAKiKl1FWWuJzRSIi+SMnw/3tZVtYd7jda8mMra8k\nEs7JUkVEclLOJebxk2eYMnWhF+y1VTHi5RFNNSAiMgA5d537g7+dz9ozFQCEQ0HG1lcx6e4bfa5K\nRCS/5NTIfc3mPSzeddwbtY8fnaSsNOd+/oiI5LycCXfnHI+/tMAL9kQsSn0ypnaMiMgg5MyweNHq\n7aw/dBYCYcxg3Oik2jEiIoOUEyP3zs4UL7yxjLZAOZB+fF68TJc+iogMVk6E+9wlm9hzsI0Uhln6\nhiW1Y0REBs/3tkzbh6d4aNoqDmRuWBpVkyAcCupuVBGRIfB15H6uvYP/99gbHGgvIYURDgUZVR3X\nDUsiIkPka4r+5Pn3mLsP7wqZCZlLH9WSEREZGt/aMq0H2nh15X4v2EfXJqhLxnSFjIhIFvRr5G5m\nE81sk5ltMbMHe9nnYTNrMbNVZvaJSx3v+Mkz/PzxGXS69HpJOMiEhiqN2EVEsqTPcDezAPAIcBtw\nPXCXmV3XY58vA1c6564G7gV+09vxzrV38PePvsaGw+cfnXf52Boe/daniuokanNzs98l5Ax9F+fp\nuzhP38XQ9GfkfgvQ4pzb5ZxrB54H7uixzx3AUwDOuSVAwszqL3aw/zlpDnP3wfuB9OPyLhtVRU2i\nbLD15y39xT1P38V5+i7O03cxNP0J9zFAa5f1PZltl9pn70X2Aeg2d0xDbZy6ZFztGBGRLBvxE6of\nBXtlRZQrxlRz+ydHF1U7RkRkJJhz7tI7mN0KNDnnJmbWfwg459zPuuzzG2Cec+6FzPom4D865w72\nONalP0xERC7KOWcD2b8/I/elwFVmNh7YD9wJ3NVjn+nAfcALmR8Gx3oG+2CKExGRwekz3J1znWZ2\nPzCLdI9+inNuo5ndm37ZTXbOvWFmXzGzrcBJ4FvDW7aIiFxKn20ZERHJPyM2/UB/boQqBmY21szm\nmtl6M1trZt/3uyY/mVnAzFaY2XS/a/GbmSXM7EUz25j5+/EnftfkBzP7X2a2zszWmNkzZlZU83+b\n2RQzO2hma7psqzKzWWa22cxmmlmir+OMSLj350aoItIB/G/n3PXAp4H7ivi7AHgA2OB3ETniV8Ab\nzrmPAR8HNvpcz4gzs9HA94AbnXM3kG4d3+lvVSPuCdJZ2dUPgdnOuWuBucDf9HWQkRq59+dGqKLg\nnDvgnFuVWT5B+h/wRe8JKHRmNhb4CvC437X4zcziwH9wzj0B4JzrcM596HNZfgkC5WYWAsqAfT7X\nM6KccwuAth6b7wCezCw/CXy1r+OMVLj350aoomNmE4BPAEv8rcQ3DwH/F9CJH7gcOGJmT2TaVJPN\nLOp3USPNObcP+AWwm/TNkMecc7P9rSon1H10BaJz7gBQ19cbNHG6T8wsBrwEPJAZwRcVM/svwMHM\nbzGW+VPMQsCNwCTn3I3AKdK/ihcVM6skPUodD4wGYmb2dX+rykl9DohGKtz3AuO6rI/NbCtKmV83\nXwL+zTn3it/1+OSzwO1mth14DviCmT3lc01+2gO0OueWZdZfIh32xeZPge3OuaPOuU7gZeAzPteU\nCw5+NF+XmY0CDvX1hpEKd+9GqMyZ7ztJ3/hUrH4HbHDO/crvQvzinPuRc26cc+4K0n8f5jrn/pvf\ndfkl8yt3q5ldk9n0JYrzRPNu4FYzKzUzI/09FN2JZS78bXY6cE9m+W6gz0HhiMwt09uNUCPx2bnG\nzD4LfANYa2YrSf969SPn3Ax/K5Mc8H3gGTMLA9spwpsBnXPvmdlLwEqgPfO/k/2tamSZ2bNAI1Bt\nZruBHwP/BLxoZt8GdgF/1edxdBOTiEjh0QlVEZECpHAXESlACncRkQKkcBcRKUAKdxGRAqRwFxEp\nQAp3EZH6TuPLAAAADUlEQVQCpHAXESlA/x/LZn+fWPd4ZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f48d8a46a50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def SampleWeibull(lam, k, n=1):\n",
    "    return np.random.weibull(k, size=n) * lam\n",
    "\n",
    "data = SampleWeibull(lam, k, 10000)\n",
    "cdf = Cdf(data)\n",
    "model = pmf.MakeCdf()\n",
    "thinkplot.Cdfs([cdf, model])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise:** Write a class called `LightBulb` that inherits from `Suite` and `Joint` and provides a `Likelihood` function that takes an observed lifespan as data and a tuple, `(lam, k)`, as a hypothesis.  It should return a likelihood proportional to the probability of the observed lifespan in a Weibull distribution with the given parameters.\n",
    "\n",
    "Test your method by creating a `LightBulb` object with an appropriate prior and update it with a random sample from a Weibull distribution.\n",
    "\n",
    "Plot the posterior distributions of `lam` and `k`.  As the sample size increases, does the posterior distribution converge on the values of `lam` and `k` used to generate the sample?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise:** Now suppose that instead of observing a lifespan, `k`, you observe a lightbulb that has operated for 1 year and is still working.  Write another version of `LightBulb` that takes data in this form and performs an update. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Now let's put it all together.  Suppose you have 15 lightbulbs installed at different times over a 10 year period.  When you observe them, some have died and some are still working.  Write a version of `LightBulb` that takes data in the form of a `(flag, x)` tuple, where:\n",
    "\n",
    "1. If `flag` is `eq`, it means that `x` is the actual lifespan of a bulb that has died.\n",
    "2. If `flag` is `gt`, it means that `x` is the current age of a bulb that is still working, so it is a lower bound on the lifespan.\n",
    "\n",
    "To help you test, I will generate some fake data.\n",
    "\n",
    "First, I'll generate a Pandas DataFrame with random start times and lifespans.  The columns are:\n",
    "\n",
    "`start`: time when the bulb was installed\n",
    "\n",
    "`lifespan`: lifespan of the bulb in years\n",
    "\n",
    "`end`: time when bulb died or will die\n",
    "\n",
    "`age_t`: age of the bulb at t=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lifespan</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>age_t</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.493720</td>\n",
       "      <td>1.661899</td>\n",
       "      <td>2.155619</td>\n",
       "      <td>8.338101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.532562</td>\n",
       "      <td>5.336770</td>\n",
       "      <td>5.869332</td>\n",
       "      <td>4.663230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5.408782</td>\n",
       "      <td>9.291016</td>\n",
       "      <td>14.699798</td>\n",
       "      <td>0.708984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.055808</td>\n",
       "      <td>6.363279</td>\n",
       "      <td>10.419087</td>\n",
       "      <td>3.636721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.264431</td>\n",
       "      <td>6.724079</td>\n",
       "      <td>9.988509</td>\n",
       "      <td>3.275921</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   lifespan     start        end     age_t\n",
       "0  0.493720  1.661899   2.155619  8.338101\n",
       "1  0.532562  5.336770   5.869332  4.663230\n",
       "2  5.408782  9.291016  14.699798  0.708984\n",
       "3  4.055808  6.363279  10.419087  3.636721\n",
       "4  3.264431  6.724079   9.988509  3.275921"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "lam = 2\n",
    "k = 1.5\n",
    "n = 15\n",
    "t_end = 10\n",
    "starts = np.random.uniform(0, t_end, n)\n",
    "lifespans = SampleWeibull(lam, k, n)\n",
    "\n",
    "df = pd.DataFrame({'start': starts, 'lifespan': lifespans})\n",
    "df['end'] = df.start + df.lifespan\n",
    "df['age_t'] = t_end - df.start\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I'll process the DataFrame to generate data in the form we want for the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('eq', 0.49372023101011586)\n",
      "('eq', 0.53256153260227868)\n",
      "('gt', 0.70898391930637672)\n",
      "('gt', 3.636720876317753)\n",
      "('eq', 3.2644306885160645)\n",
      "('eq', 1.0391519214771501)\n",
      "('eq', 0.54538148700181943)\n",
      "('eq', 1.1052773283533255)\n",
      "('eq', 1.5679787046956348)\n",
      "('eq', 0.26546578499245693)\n",
      "('eq', 0.79189654373581286)\n",
      "('eq', 5.8699951990938581)\n",
      "('eq', 2.7085311291070671)\n",
      "('eq', 0.46889988780996417)\n",
      "('eq', 0.49196885895901804)\n"
     ]
    }
   ],
   "source": [
    "data = []\n",
    "for i, row in df.iterrows():\n",
    "    if row.end < t_end:\n",
    "        data.append(('eq', row.lifespan))\n",
    "    else:\n",
    "        data.append(('gt', row.age_t))\n",
    "        \n",
    "for pair in data:\n",
    "    print(pair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise:** Suppose you install a light bulb and then you don't check on it for a year, but when you come back, you find that it has burned out.  Extend `LightBulb` to handle this kind of data, too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Prediction\n",
    "\n",
    "**Exercise:** Suppose we know that, for a particular kind of lightbulb in a particular location, the distribution of lifespans is well modeled by a Weibull distribution with `lam=2` and `k=1.5`.  If we install `n=100` lightbulbs and come back one year later, what is the distribution of `c`, the number of lightbulbs that have burned out?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Now suppose that `lam` and `k` are not known precisely, but we have a `LightBulb` object that represents the joint posterior distribution of the parameters after seeing some data.  Compute the posterior predictive distribution for `c`, the number of bulbs burned out after one year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
